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Abstract 

As a part of CFD code validation efforts within the Icing 
Branch of NASA Glenn Research Center, computations were 
performed for natural laminar flow (NLF) airfoil, NLF-0414, 
with 6 and 22.5 minute ice accretions. Both 3-D ice castings 
and 2-D machine-generated ice shapes were used in wind 
tunnel tests to study the effects of natural ice as well as simu- 
lated ice. They were mounted in the test section of the Low 
Turbulence Pressure Tunnel (LTPT) at NASA Langley that 
the 2-dimensionality of the flow can be maintained. Aerody- 
namic properties predicted by computations were compared 
to data obtained through the experiment by the authors at the 
LTPT. Computations were performed only in 2-D and in the 
case of 3-D ice, the digitized ice shape obtained at one span- 
wise location was used. The comparisons were mainly con- 
centrated on the lift characteristics over Reynolds numbers 
ranging from 3 to 10 million and Mach numbers ranging 
from 0.12 to 0.29. WIND code computations indicated that 
the predicted stall angles were in agreement with experiment 
within one or two degrees. The maximum lift values 
obtained by computations were in good agreement with 
those of the experiment for the 6 minute ice shapes and the 
22.5 minute 3-D ice, but were somewhat lower in the case of 
the 22.5 minute 2-D ice. In general, the Reynolds number 
variation did not cause much change in the lift values while 
the variation of Mach number showed more change in the 
lift. The Spalart-Allmaras (S-A) turbulence model was the 
best performing model for the airfoil with the 22.5 minute 
ice and the Shear Stress Turbulence (SST) turbulence model 
was the best for the airfoil wuth the 6 minute ice and also for 
the clean airfoil. The pressure distribution on the surface of 
the iced airfoil showed good agreement for the 6 minute ice. 
However, relatively poor agreement of the pressure distribu- 
tion on the upper surface aft of the leading edge horn for the 
22.5 minute ice suggests that improvements are needed in 
the grid or turbulence models. 

Nomenclature 

LTPT Low- Turbulence Pressure Tunnel (NASA Langley) 


MVD Median Volume Diameter (in \sm ) 

LWC Liquid Water Content (in g/in ) 

IRT Icing Research Tunnel (at NASA Glenn) 

AOA Angle-of- Attack 

S-A Spalart-Allmaras Turbulence Model 

SST Shear Stress Transport Turbulence Model 

C p Pressure coefficient 

C, Maximum 2D lift coefficient 

1 max 

a C{ Angle where maximum lift occurs 

Introduction 

NLF airfoils were designed to achieve lower cruise drag 
coefficients while retaining the high maximum lift coeffi- 
cients intended for use in general aviation. Experimental 

worlJ^ using a clean wing was performed at the NASA Lan- 
gley LTPT tunnel in the 1980s but the aerodynamic measure- 
ment of icing effects has not been performed for this type of 
airfoil even though a minor ice accretion on such a wing 
could result in a significant penalty in the aerodynamic per- 
formance. In the aerodynamic measurement of iced airfoils 
and wings performed in the past, many researchers have used 
ice models w ith either very smooth surfaces, ones with sand 
grain roughness or bead type protrusions to simulate the real 
iced surface. These shapes were assumed to provide perfor- 
mance degradation reasonably close to real ice but the actual 
difference has not been verified. 

The primary purpose of this research was to evaluate the 
performance of the CFD code against the experimental data. 
An additional purpose of the test and computation was to 
investigate the effects of Reynolds numbers as well as Mach 
numbers on the performance of the iced airfoil. In order to 
perform this study, actual cast ice made from molds of ice 
formed in the Icing Research Tunnel (IRT) at NASA Glenn 
as well as machine generated two-dimensional projections of 
the ice shapes obtained at one spanwise location w^ere used 
to investigate the effects of smoothing on the aerodynamic 
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characteristics. Comparisons were also made for the lift 
characteristics including the maximum lift and stall angles 
which are very important factors for safe operations of air- 
craft. 

Approach 

The ice shapes used in this study were accreted on the NLF- 
0414 airfoil in the IRT under the following icing conditions. 
Icing spray conditions which produced the ice shapes were a 
median volume diameter (MVD) of 20pm , a static air tem- 
perature of -5°C , a model attitude of 2° , a liquid water con- 
tent (LWC) of 0.54 g/m 2 3 and 6 and 22.5 minute spray times. 
Castings were made of these ice shapes for the LTPT aero- 
dynamic test. The Mach number was 0.21 and Reynolds 
number was 4.6 million based on the chord. The numerical 
grid sensitivity analysis explained later in this study was 
based on this experimental condition. 

Wind Tunnel TesP ^ 

The NLF-0414 model tested in the LTPT had a chord of 36 
inches and a span of 36 inches. It had a removable leading 
edge in order to attach different ice shapes. Figure 1 shows 
the three-dimensional cast ice pieces mounted with a stain- 
less steel piece for pressure taps sandwiched between them. 
Figure 2 shows a similar mounting of the machine-generated 
two-dimensional model. Inside the test section of the LTPT, 
side wall boundary layer control was used to ensure the two- 
dimensionality of the flow'. The experimental lift data was 
obtained at all angles including the post-stall angles but due 
to the enormous time requirements to operate the wake sur- 
vey, drag data was obtained only at limited number of angles 
for most cases. Both force balance and integration of surface 
pressure data were used for the calculations of global force 
characteristics. Due to its accuracy, the surface pressure 
was used to obtain the overall lift and the force balance data 
showed good agreement with it. Full experimental results 
can be found in reference 2. Low speed wind tunnel wall 
and blockage corrections^ were applied to the experimental 
data and these were compared to computational results run 
under external flow' conditions. Advantages of LTPT were 
wider range of Reynolds and Mach numbers and low turbu- 
lence level than ones available in smaller atmospheric tun- 
nels. The test Reynolds numbers were ranged from 3 to 10 
million and the Mach numbers were ranged from 0.21 to 
0.29. 

Sutface Modeling 

Past research efforts showed that proper modeling of the iced 
surface and a grid sensitivity study are essential for the accu- 
rate prediction of aerodynamic properties^. When the sur- 


face of the 3-D ice castings was modelled for 2-D 
computations, the digitized ice shape at one spanwise loca- 
tion was used. Only the areas which made grid generation 
almost impossible were treated to have smoother geometry 
fit for grid generation. Surface modeling tool kit ‘Srnag- 
glce^’ developed in-house at NASA Glenn was used for the 
treatment of such areas and to construct smoothed geome- 
try for the fabrication of two-dimensional ice shapes by Ste- 
reo-lithography machine. The same approach reported in 
references 4 was used to model the iced surfaces by using 
20 7c level of control points for the 22.5 minute ice and 50 7c 

for the 6 minute ice^. 

Grid Generation 

The modeling of complicated 3-D surfaces is still technically 
very challenging and the computation of such cases would 
be very expensive. Although the ice castings represent truly 
three-dimensional surfaces, all grids for this study were gen- 
erated for 2-D computations. Figures 3a and 3b compare the 
original digitized ice shape and the smoothed ice shape for 
the 22.5 minute ice. A similar smoothing was also performed 
for the 6 minute ice. A series of grids with variable densities 
both in circumferential and normal directions were first gen- 
erated to determine the optimum resolution. Figures 4a 
shows the grid block strategy around the leading edge of the 
22.5 minute ice. A C-type grid was used for the inner block 
with the farfield boundaries placed at a distance of 15 chord 
lengths from the body surface in all directions. This block- 
ing method was used to ensure good grid quality control for 
the inner block around the complicated ice shapes. Figure 
4b shows grid lines near the leading edge of the smaller 6 
minute ice. Elliptic grid generation was performed using 

the GRIDGEN program. 

Numerical Methods 

The WIND code^ capable of simulating complex flow 
fields was used for the calculation of the current problems. It 
is capable of multi-zone, 2-D/3-D, implicit/explicit, and 
steady/unsteady computations with various turbulence mod- 
els. The governing equations in WIND code are formulated 
using the finite-volume approach and the specification of the 
discretization on the right hand side is flexible. The choice 
includes standard Roe upwind, physical Roe upwind for 
stretched grids and Coakley upwind. The spatial accuracy 
can be adjusted to the fifth order if desired by the user. Sev- 
eral turbulence models including the Spalart-Allmaras (S- 

A) [7] and the Baldwin-Barth (B-B) f8] one-equation turbu- 
lence models and Shear Stress Transport (SST) two-equation 

turbulence model^ were used for the grid sensitivity test for 
the determination of the best performance. The SST model 
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is a hybrid model which blends the solution of the k - 03 
model near a solid surface and the k - £ model elsewhere 
including the shear layer. A fully turbulent flow was assumed 
for the computation considering the roughness of the ice 
shapes and partly due to lack of ability to include transition 
effects. At the far field boundary, a non-reflecting type 
boundary condition was applied. 

Grid Sensitivity' Test 

In order to achieve high quality numerical results, a series of 
grids for the inner block having different densities in both 
normal and circumferential directions were constructed (the 
dimension of the outer block was fixed). The WIND code 
was then run with three different turbulence models 
described above, using these grids, and the aerodynamic 
properties were compared to the experiment. 

The first set of grids constructed were used for an investiga- 
tion of the normal direction sensitivity, which was followed 
by a study of the effect of packing grid points in the circum- 
ferential direction. Combination with various turbulence 
models was also attempted to determine the proper model for 
each ice shape. In the first study applied to the 22.5 minute 
3-D ice, the number of grid points used in the normal direc- 
tion were 71, 81, 91, and 101 while the circumferential 
direction grid number was fixed at 491. Figure 5 shows the 
effect of this normal direction grid packing. It showed that 
about 101 grid points were needed for good agreement with 
experimental data When the number of grid points in normal 
direction was further increased from 101, it did not result in 
significant change of lift. Application of the three turbulence 
models for 491 x 81 inner block grid showed that the Spal- 
art-Allmaras models was the best one for this type of rela- 
tively large ice with horns (Figure 6). 

Based on this result, the normal direction grid was fixed at 
101 points and the circumferential grid densities of 321, 361, 
421, 441, 491, and 501 were applied. Figure 7 shows the 
result with three grid resolutions indicating that the change 
of grid density in circumferential direction did not affect the 
lift as much as the normal direction change did. Also in this 
case, grid points more than 441 did not change the maximum 
lift value significantly. In all cases, the stall angle was the 
same regardless of the number of points used. The reason- 
able resolution of 441x101 was used for the rest of the com- 
putation for the 22.5 minute 3-D ice. When the same 
approach was applied to the 22.5 minute 2-D ice, 401x101 
grid was determined to be the best one. In the case of 6 
minute ice, the best performing turbulence model was SST 
model and the grid resolution for the 2-D and 3-D ice w 7 ere 
351 x 81 and 371 x 81 respectively. For all these studies, the 
non-dimensional first grid spacing in the normal direction 


(yl, or minimum wall spacing) was 4.0 x 10 6 which was 
determined from numerical experiments. This value resulted 
in average value of y+ less than 1 .0 for all ice shapes and 
flow conditions. 

Discussions of Results 

Total of five different flow conditions were selected for CFD 
runs for 6 minute 2-D and 3-D ice (table 1) while seven con- 
ditions were applied for the 22.5 minute ice (table 2). Out of 
these flow conditions, Mach number of 0.21 and Reynolds 
number of 6.4 million was selected as the representative con- 
dition for the discussions to follow. 

6 minute 3-dimensiona I cast ice 

Figure 8a shows a comparison of lift versus AOA curves 
between the computation and the experiment for Mach num- 
ber of 0.21 and Reynolds number of 6.4 million. The maxi- 
mum lift obtained by experiment w-as 1.023 and the 
computational prediction w 7 as 1.019, which w 7 as a difference 
of 0.49F But the computationally predicted stall angle was 
8.3 degrees which was one degree lower than that of the 
experiment. The difference in the maximum lift ranged from 
0.4 percent for this condition to as much as 6.5 percent for 
M=0. 12 and Re = 6.4 million. The stall angle remained 
unchanged at 8.3 degrees for all five flow 7 conditions and the 
experimental stall angles were all at 9.3 degrees. Relatively 
good agreement was also observed in the comparison of the 
C for all the cases. At the same flow condition as Figure 

8a, comparisons at the angle of attack of 0.1 degree and 8.3 
degrees (Figures 8b and 8c) showed generally good agree- 
ment except near the trailing edge. At 8.3 degree AOA, the 
C by computation showed slight deviation from the experi- 
ment on the upper surface. It should be noted that some dif- 
ferences might come from the fact that the experimental 
pressure is an uncorrected one obtained at 8.0 degrees (pres- 
sure coefficients cannot be corrected while the global forces 
can be) but the computation was performed at the corrected 
angle of 8.3 degrees for the fair comparison of lift values. 
Considering this level of agreement of C p curves, the com- 
putationally predicted separation region on the upper surface 
and the velocity distribution pictured in Figure 8d is a rea- 
sonable representation of the experimental condition though 
velocity was not measured in the tests. The next two figures 
9 and 10 summarize the effect of Reynolds number and 
Mach number on the lift. The fact that the Reynolds number 
did not have as much influence as the Mach number w 7 as in 
agreement with experiment but the margin of difference was 
greater than experiment especially for the Mach number 
effect. The decrease in C u w 7 as 1.59c in going from Mach 
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number of 0.12 to 0.21 and 8.8% when going from 0.12 to 
0.29. In comparison, the experiment showed 1.5% and 
3.2% changes. This suggests that the physical model in the 
code may be sensitive to the compressibility effect caused by 
the increase in Mach number. 

6 minute 2-dimensional smooth Lee 

For the machine-fabricated smooth ice, comparison showed 
somewhat larger variations than the 6 minute 3-D ice case. 
At Mach number 0.12, the a C) for experiment was 10.3 

degrees while the computational result indicated that it was 
9.3 degrees. For the flow 7 condition of Mach number of 0.21 
and Reynolds number of 6.4 million, the largest difference 
was observed between the experiment and the computation 
(Figure 11). The a c< for the computation occurred 2 

degrees lower and C /mw was 2.1% lower but for all other 
flow 7 conditions, one degree difference in the a C/ w 7 as 

observed as in the 3-D case. The overall difference in C Ut 
for the five conditions increased only slightly compared to 
the 3-D cast ice case. Similar trend in the influence of Mach 
number and Reynolds number on the lift w 7 as also observed. 
The C imv was 9.5% higher for the 2-D ice shape than for the 
3-D ice shape as shown in Figure 12a. This is in good agree- 
ment with experiment since its margin was 9%J^ Figure 
12b compares the C between the 2-D and the 3-D cases and 
it shows some differences on the upper surface. In general 
for both 2-D and 3-D cases, computations predicted a C( 

and C t somew'hat conservatively by producing slightly 
lower values. 

22.5 minute 3 -dimensional cast ice 

Figures 13a through 13g show 7 lift, C p , and flow field char- 
acteristics for the larger 22.5 minute 3-D cast ice under the 
representative flow condition. In this case, computationally 
predicted a C/ w r as one degree higher at 6.2 degrees than the 

experimental result of 5.2 degrees. Computational C tmax was 
2.0% lower than that of experiment. Although the lift pre- 
dicted by the computation showed relatively good agreement 
with the experiment, the C curves compared at 0. 1 degree 

and 6.2 degree AOAs indicated that the computational pre- 
diction did not match the experiment at both angles. Espe- 
cially at AOA of 6.2 degrees, the computational C p showed 

much faster pressure recovery 7 on the upper surface than the 
experiment. This would cause differences in the separation 

bubble size and the velocity distribution in it. Bragg, et al™ 
performed experiments to determine boundary layer transi- 


tion for a NACA 0012 airfoil with simulated icing scale 
roughness at low Mach numbers and Reynolds numbers 
ranging from 0.75 to 2.25 million. They showed that the 
fully developed turbulent profiles were not measured until 
approximately 40% chord location. Although neither flow 
conditions nor roughness elements on the surface of the 
study match those of the current study, it is enough to sug- 
gest that a certain level of transition effect should be intro- 
duced in the CFD computations. Apart from these pressure 
distribution differences, the Reynolds number and Mach 
number variations showed very similar effect as the 6 minute 
ice cases. Figures 1 4a and 14b showed that increasing the 
Reynolds number did not have significant effect on the lift at 
Mach numbers of 0.12 and 0.21 while Figure 1 4c show r ed 
that the Mach number change had more significant effect. 
The decrease in C /wj was 9.6% in going from Mach number 
of 0.12 to 0.21 and 12.3% when going from 0.12 to 0.29. 
This is a quite a bit of contrast to the experiment which 
showed only about 0.5% change in the C tmax for the three 
Mach numbers tested at Reynolds number of 6.4 million. 

22.5 minute 2-dimensional smooth ice 

For this ice shape, the largest lift differential was observed at 
a Mach number of 0.21 and Reynolds number of 6.4 million 
as shown in Figure 15a. The a c . by experiment was one 

degree higher at 7,3 degrees than the computation showing 
one degree shift of the stall angle compared to the 3-D case. 
Figure 15b showed similar differences in C n as the 3-D case 

but the pressure difference aft of the upper horn and the 
lower surface was less than that of the 3-D cast ice case. Fig- 
ure 16a and 16b indicated that there was 2.8% decrease in 
the Ci max going from 2-D to 3-D and the C p data showed no 

significant difference. This is the greatest disagreement 
among the 7 conditions (see table 2) computed but at the 
lower Mach number of 0.12, better agreement with the 
experiment was observed as shown in Figures 17 and 18. At 
Reynolds number of 3.0 million, the a C; matched at 7.2 

degrees and the computational C lmar was only 0.95% lower 

than the experiment. At higher Reynolds number of 6.4 mil- 
lion, the stall angle was one degree higher for the experiment 
(8.2 degrees versus 7.3 degrees) and the computational C lmM 
was 2.4% lower than the experiment. 

Effect o f 2-D and 3-D ice on the lift 

Figures 19 and 20 summarize the effect of 2-D and 3-D ice 
on the lift. Compared to the computational lift for the clean 
airfoil, accretion of 3-D ice resulted in 37% drop of lift for 6 
minute ice and 50% drop for 22.5 minute ice. Similarly for 
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the 2-D ice, 31 7c and 487 drop of lift was observed respec- 
tively. 

Concluding Remarks 

A parametric computational study was performed to evaluate 
the ability of WIND code to predict the performance degra- 
dation of NLF-0414 wing with realistic 3-D cast ice and 2-D 
machine-generated ice shapes. Computations were per- 
formed in 2-D domain with digitized geometry for the 3-D 
ice shape and smoothed geometry for the 2-D ice shape for a 
range of Mach and Reynolds numbers. The WIND code gen- 
erally predicted conservatively lower lift and the a Ci were 

in agreement with the experiment within one or two degrees. 
The r values were in good agreement for the 6 minute ice 

shapes and the 22.5 minute 3-D cast ice, but were somewhat 
lower in the case of 22.5 minute 2-D ice for Mach numbers 
greater than 0.21. Both computations and experiments indi- 
cated that Mach number had more effect on the lift than Rey- 
nolds number did although greater change on the lift was 
observed for computations with the increased Mach number. 
The Spalart-Allmaras turbulence model was the best model 
for the 22.5 minute ice and the SST model was the best for 
the 6 minute ice and the clean airfoil. The pressure distribu- 
tion on the surface of the iced airfoil showed good agreement 
for the 6 minute ice but relatively poor agreement on the 
upper surface aft of the leading edge horn for the 22.5 
minute ice. 

The following observations should be noted for future stud- 
ies similar to this one. First, agreement in lift characteristics 
(or global forces) does not guarantee the agreement in pres- 
sure distributions. Second, more comprehensive experimen- 
tal data including drag, pitching moment, and velocity 
measurements would be helpful for an accurate assessment 
of the capabilities of the CFD code. The grid sensitivity 
study and the selection of the proper turbulence models thus 
could be performed mainly by comparing the lift values 
against the experiments limiting the accuracy of such a 
study. Third, the differences in the pressure distribution on 
the surface of the airfoil suggest that more investigations are 
needed in constructing better grids especially in the shear- 
layer-dominated regions aft of large horns, grid adaptations, 
and better turbulence models which can effectively account 
for the existence of boundary layer transition. Fourth, con- 
sidering the fact that the SST model was the best for the 
smaller 6 minute ice shape but S-A model was the best for 
the larger 22.5 minute ice shape, more study is recom- 
mended for the determination of the best performing turbu- 
lence models for different sizes and shapes of ice. Fifth, 
further investigation is needed to explain the exact causes for 
the differences in the global force characteristics and the 


pressure distributions between the computations and experi- 
ments as shown in the current study. 
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Table 1. Flow conditions used for the computations of the 6 
minute ice 


Re xlO' 6 

M - 0.12 

M = 0.21 

M=0.29 

4.6 


X 


6.4 

X 

X 

X 

10.0 


X 



Table 2. Flow conditions used for the computations of the 
22.5 minute ice 


Re x 10 ' 6 

3.0 
4.6 
6.4 

10.0 


A/= 0.12 


X 


M = 0.29 



Figure 1 3-D ice castings mounted in the test section 




Figure 2 2-D smooth ice mounted in the test section 


Figure 4b Grid around the leading edge of the 6 min. ice 
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-igure 8b Comparison of Cp for 6 min. 3-D ice ; AOA = 0.1 deg 


Figure 8d Flow pattern around the ieadina edge of 
the 6 minute ice at AOA = 8.3 degrees 


M = 0.21, Re = 6.4 E-t-6, WIND code with SST turb. model 

-6 r « 1 ■ 1 ' 



1.4 


M - 0.21 , WIND code with SST turb. model 


AOA (deg.) 


Figure 8C Comparison of Cp for 6 min. 3-D ice : AOA = 8.3 deg. 


Figure 9 Reynolds number effect on the 6 min. 3-D ice 
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Re = 6.4 mil.. WIND code with SST turb. model 



Figure 1 0 Mach number effect on the 6 min. 3-D ice 


M = 0.21 . Re * 6.4 E+6 Wind code with SST lurb. model 



Figure 1 1 Lift Comparison for the 6 min. 2-D ice shape 


M = 0 21 , Re = 6 4 E+6, Wind code with SST turb. model 



Figure 1 2a Lift Comparison for the 6 min. 2-D & 3-D ice 


6 min. ice, M = 0.21 , Re = 6.4 E+6, WIND code with SST 



Figure 12b Cp comparison for 2-D and 3-D ice : AOA = 8.3 deg. 
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M = 0.12, Re = 3.0 E+6. Wind code with S-A lurt) . model 


3-D rough ice : M = 0.21, Re = 6.4 mil., WIND code 




: igure 1 7 Lift Comparison for the 22.5 min. 2-D ice shape Figure 1 9 Effect of the 3-D ice on the lift 


M = 012 Re s 6.4 E+6 Wind code with S-Aturb. model 2-D smooth ice : M = 0.21, Re = 6.4 mii , WiND code 




Figure 1 8 Lift Comparison for the 22,5 min. 2-D ice shape 


Figure 20 Effect of the 2-D ice on the lift 
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